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Abstract 

Swimming micro-organisms such as bacteria or spermatozoa are typically found in dense 
suspensions, and exhibit collective modes of locomotion qualitatively different from that dis- 
played by isolated cells. In the dilute limit where fluid-mediated interactions can be treated 
rigorously, the long-time hydrodynamics of a collection of cells result from interactions with 
many other cells, and as such typically eludes an analytical approach. Here we consider the 
only case where such problem can be treated rigorously analytically, namely when the cells 
have spatially confined trajectories, such as the spermatozoa of some marine invertebrates. We 
consider two spherical cells swimming, when isolated, with arbitrary circular trajectories, and 
derive the long-time kinematics of their relative locomotion. We show that in the dilute limit 
where the cells are much further away than their size, and the size of their circular motion, 
a separation of time scale occurs between a fast (intrinsic) swimming time, and a slow time 
where hydrodynamic interactions lead to change in the relative position and orientation of the 
swimmers. We perform a multiple-scale analysis and derive the effective dynamical system — of 
dimension two — describing the long-time behavior of the pair of cells. We show that the system 
displays one type of equilibrium, and two types of rotational equilibrium, all of which arc found 
to be unstable. A detailed mathematical analysis of the dynamical systems further allows us 
to show that only two cell-cell behaviors are possible in the limit of t — > oo, cither the cells are 
attracted to each other (possibly monotonically), or they are repelled (possibly monotonically 
as well), which we confirm with numerical computations. Our analysis shows therefore that, 
even in the dilute limit, hydrodynamic interactions lead to new modes of cell-cell locomotion. 

Keywords: Hydrodynamic interactions - Swimming cells - Collective locomotion - Multiple- 
scale analysis 
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1 Introduction 



Micro-organisms such as bacteria and simple eukaryotes are found in nature in a variety of environ- 
ments, from large water masses (ocean, lakes, rivers) to the fluid components of plants and animals. 
In all, they represent half of the world's biomass, and have therefore major biological consequences 
on the health and survival of most other organisms. 

When a micro-organism has the ability to swim in a viscous fluid, then its motion is the 
complicated result of the local transport by the moving fluid it resides in, and of its intrinsic 
swimming. Given the small size, £, of these micro-organisms (typically i ~ 1-10 fim) and the 
small swimming velocities, V (typically V ~ 10 — 100 fim/s), the Reynolds number, Re = V£/v, 
is much smaller than 1 (here v is the kinematic viscosity of the fluid). For such swimmers, the 
interactions with the surround i ng flu id are therefore dominated by viscous stresses, and inertial 



effects are negligible (jLighthill 



19751). As a results, the velocity and pressure fields around the 



swimmer satisfy Stokes' equations (jHappel and Brenner . 



1965; 



Kim and Karillal . 



1991 



Most classical wor 
individual organisms ( 


<; on the dynamics of swimming cells considered the mechanics and phvsics of 


Lie 


;hthill 


1976; 


Brennen and Winet. 


1977; 


Blum and Hines. 


1979; 


Childress. 


1981; 


Lauga and Powers . 


20oJ; 


Brav. 


2000 


). However, cells are typically found in large dense sus- 



pensions, and display collective modes of locomotion which are qualitatively different from that of 
individual cells. For example spermatozoa populations can be as large as millions, and in some 
spec ies display aggregat ion and cooperative l ocomotion. Such is the cas e for wood mouse spermato - 



zoa (Moore et al. 



20021 ) , as well as opossum (IMoore and Taggart 



19951 ) and fishfly (IHayashi 



19981 ) 



Concentrated bacterial suspensions display large-scale coherent and inter mittent collective swim- 



ming, with length and velocity scales much larger than that o f a sin gle cell (jMendelson et al 



Dombrowski et al 



2004 



Sokolov et al 



diffusion of suspended particles (jWu and Libchaber 



2007 



Cisneros et al 



2000 



1999 



20071') . and re s ulting in an enhanced 



Kim and Breuerl . 



20041 ). 



Significant work has been devoted to the theoretical modeling of collective effects in cell locomo- 



tion. Building on early work s 



rowing that dipo 



ming cells lead to aggregation (jGuell et al. 



e-dipole hydrodynamic interactions between swim- 



19881 ). two distinct approaches have been considered. On 



one hand, continuum studies have been proposed in the dilute limit. Classical work on bioconvec- 
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tion neglected the pres e nce o f swimming cells altogether ({Childress et al. 



1992 



Hill and Pedley 



1975 



Pedlev and Kessler 



20051 ). When the swimmer size is small compared to the typical inter- 



swimmer distance, the first effect of a self-propelled micro-organism is to modify t he local stresses i n 



the flow by creating a local dipolar (or stresslet) forcing on the surrounding fluid (jBatchelor 



1970). 



Within this framework, studies have discove red long-wavelength hyd r odynamic instabilities occur 



ring 



in suspensions of self-propelled bodies (jSimha and Ramaswamvl . 



2002 



Saintillan and Shelley 



20081 ). The resulting nonlinear state, some times referred to as "bacterial tur 



3ulence" has also been 



reproduced using continuum simulations (jAranson et al. 



2007 



Wolgemuth 



20081 ). On the other 



hand, a number of studies have focused on the discrete nature of the "N-swimming body" problem, 
and solved numerically for the dynamics of each self -propeled body. Models of increasing complex : 



ity h ave represented the swimmer as a point- dipole 



20081 ). a line distribution o f surface stress 



tion of tangential velocity (jlshikawa and Pedley 



Hernandez-Ortiz et al. 



Saintillan and Shel 



20071) 



a 



2005; 



Underhill et al. 



20071 ). or a surface distribu- 



20081 ) . and have reproduced some of 



the instabilities, diffusive behavior, and nonlinear dynamics observed experimentally (see also 



Mehandia and Nottl . 



20081 ) . The subtle role of hydrodynamic interaction s in al 



of lo comotion was also recently pointed out ([Alexander and Yeomans 



2008 



owing for new modes 



Lauga and Bartold . 



2008). In parallel, work in the physics community has discovered phase-transitions to collective 



motion in kinematics mod els of large populat i ons of self-propelled bodies without the need for 



hydrodynamic interactions (IVicsek et al 



1995; 



Czirok et al 



1997 



Gregoire and Chate 



200J) 



From a theoretical standpoint, collective locomotion is a difficult problem. To be treated satis- 
factorily, the motion of N 3> 1 identical swimmers should be integrated in time. In the dense limit, 
no simple model is available to correctly describe the interplay between hydrodynamic and steric 
(excluded-volume) interactions. One simplification is to consider the dilute limit, in which hydro- 
dynamic interactions can be described by dipole-dipole interactions, and steric interactions can be 
neglected. However in this limit, hydrodynamic interactions are weak, and an order-one change 
in the trajectory of a straight-swimming body can only result from a large number of successive 
interactions with different swimmers. In other words, even in the dilute limit, one needs in general 
to study N ^> 1 cells to quantitatively capture their coupled dynamics. 

In this paper, we consider the only situation in which the case of N = 2 swimmers can give rise 
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to order-one changes in the long-time limit of their positions and orientations even in the dilute 
limit, namely when the individual swimmers have spatially confined intrinsic trajectories. In that 
case, even small hydrodynamic interaction can accumulate over times long compared to an intrinsic 
swimming time, and lead to nontrivial nonlinear dynamics of the coupled system. By studying in 
the long-time limit one of these prototypical situations, we hope to obtain important physical and 
mathematical insight on the general behavior of larger populations. 

We focus our study on the particular situation where the intrinsic motio n of the micro - organ isms 



is circular. This is the ca se, for example, for sea urchin spermatozoa (jRiedel et al. 



other marine invertebrates ([Goldstein! 



2003 ). or 



19771 ). We consider two identical but arbitrary model cells, 
and assume they are widely separated. This assumption allows us to propose a simple general 
representation of cell-cell hydrodynamic interactions in £j2j We then show that a separation of 
time scales occurs, with a short time representing the intrinsic swimming time for each cell, and 
the long time being the time one has to wait for repeated hydrodynamic interactions to lead to 
order-one changes in the swimmers trajectories. This separation of time scales allows us to perform 
a multiple-scale analysis of the coupled dynamics in The equilibrium configurations of the two 
cells, as well as their stability, are studied in §U The time-averaged equations are reduced to a 
two-dimensional dynamical system whose behavior is analyzed in detail. In particular, we show 
that only two long-time behaviors can arise, as determined solely by the initial relative orientations 
of the swimmers: Either hydrodynamic interactions have a net repulsive effect and the swimmers 
eventually swim infinitely far away from each other, or they have a net attractive effect, and lead 
to collisions (or aggregation) of the two swimmers. Any relative equilibrium or limit cycle is found 
to be unstable, and we therefore do not observe any organization of the swimmers' motion through 
hydrodynamic interactions. Our modeling assumptions, some possible extensions, and the relevance 
to biological locomotion are discussed in £}5j 
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2 Equations of motion of two Stokesian swimmers 
2.1 Intrinsic motion 

We first consider an isolated swimmer, whose intrinsic motion is the superposition of a translation, 
Uoe, and a rotation, Qo e ' i where e and e' are two directions rigidly attached to the swimmer. We 
neglect here the shape changes of the swimmer, assuming the swimming motion is generated by 
surface displacements that are small compared to the general dimensions of the swimmer. This is 
the case for example for a so-called sq uirmer with tangentia l displacements for which the s hape is 



at all time a sphere of constant radius (llshikawa et al 



2006; 



Ishikawa and Pedley 



20071) 



a 



20081) 



The two directions e and e' are fixed in the frame attached to the swimmer and their relative 
orientation is independent of time. In the absence of Brownian motion, the resulting equations of 
motion for the model cell are given by 

dr de „ , de' . , 

- = E7„e, - = O e<xe, ^ = (1) 

Considering only the non-trivial case where Uq 7^ 0, three situations can be considered: 

— Slo = 0: If isolated, the swimmer keeps a fixed orientation and swims along a straight line at 
constant speed 

— f?o 7^ and e • e' = 0: The isolated swimmer has a periodic motion along a circular trajectory 
of radius Uo/£Iq normal to e' and the period of the motion is 2ir/Qo. 

— General case: When 7^ and e • e' 7^ 0, the swimmer trajectory is an helix (right-handed 
if e • e' > 0, left-handed otherwise). The pitch of the helix is {2ttUq/VLq) e • e', the radius of 
the circular projection is Uq/SIq\/1 — (e • e') 2 

In this paper, we consider the case of swimmers with circular trajectory, so that e • e' = (see 
illustration in Fig. [T]). 
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Figure 1: Isolated rotating swimmer in a circular trajectory. The intrinsic velocity of the swimmer 
is the superposition of a translation parallel to e and a rotation along e'. Here, it is assumed that 
e • e' = which leads to a circular trajectory of radius p = Uq/VIq. The local basis (e, e',e x e') 
moves rigidly with the swimmer. 



2.2 Far-field velocity and vorticity field created by a general swimmer 



In this work, we propose a study of hydrodynamic interactions in the far-field limit, considering 
only the dominant contribution to the velocity field setup by the self-propelled bodies. The advan- 
tage of such an approach is to avoid having to focus on one particular geometry and gait of the 
swimmer considered. More detailed studies of hydrodynamic interactions ca n be obtained by con 
sidering the full flow field created by a bio l ogical ly realistic self-propelled cell (jlshikawa et al 



Ishikawa and Hotal . 



2007 



Gyrya et al 



2006; 



Ishikawa et al 



20071 ) or for simplified swimmer models (jPoolev et al 



2006 



2009). In the former case, the flow field must in general be solved for nu- 



merically, while in the latter, the simplification of the geometry and swimming stroke allows for 
analytical treatment. 

In general, a self-propelled cell creates its intrinsic swimming velocity, U, and angular velocity, 
CI, by imposing a displacement of its surface. This is the case for all well-studied motile cells, 
including spermatozoa, bacteria, ciliates and algae. We denote the swimmer surface S fi . This 
stroke velocity field is noted u (s) with s the position vector as measured from a point within or in 
the vicinity of the body position, and fixed in the absolute reference frame. The absolute velocity 



6 



at the boundary of the swimmer can therefore be written 



u(s)=U + Oxs + u y (s), forse^. (2) 

Let u(x) be the velocity field resulting from this swimming pattern and cr the associated stress 
field, so that f = cr ■ n is the force per unit area applied by the fluid on the swimmer's boundary 
with n the normal unit vector pointing into the fluid domain. The fluid velocity field u at a point x 



outsi de the swimmer can be expressed using the single-layer and double-layer potentials (jPozrikidia . 



19971 ) 



~[ /,(s)G^(x,s)dS(s)--L / «f(s)T iifc (x,s)n fc (s)dS(s), (3) 

SlTfl Jy »7T Jy 

where Gy(x, s) is the Green's function corresponding to the flow field at x generated by a unit 
point force located in s, and l£,-&(x, s) is the corresponding stress tensor, and where we have used 
Einstein's summation notation in Eq. ([3|). The tensors Gij and Tij k are the Green's function and 
corresponding stress tensor for the free flow case, 

Gij ( X ,s) = S f + r -^>, T yfc (s,x) = -^p, withr = x-s. (4) 

In the far-field approximation, |x| ^> |s|, and expanding Eq. ([3]) in Taylor series, the flow field is 
obtained as 

= .G^O) I' /i(a)d5(s) _ J_ ^( X)0 ) f fr(s) Sk dS(s) (5) 
T ijk (x,0) f u f {s)nk{s)dS{s)+0 ( a3 



8vr jy 

with a the typical size of the swimmer. From Eq. (jj]), we get 

dGjj _ dGjj _ SjfTk - S ik rj - 5 jk rj ^ 3rjrjr k 
ds k dr k r 3 r 5 

When Re = 0, the inertia of the swimmer is negligible and the total force and torque applied by 
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the fluid on the swimmer must vanish; therefore 



/i(s)dS(s) = 0, 



(7) 



y 



and jy, fi(s)skdS(s) must be a symmetric tensor. Consequently, the first term in Eq. © vanishes, 
and only the symmetric part in i and k of dGij/ds^, obtained in Eq. (J6j), must be retained in 
Eq. ©. The second and third terms in Eq. ([5]) behave like 1/r 2 far from the swimmer: The 



dominant velocity field f ar from the swimming body is dipolar, and dominated by a so-called 

~3) 



stresslet ( Batchelor 



197C 



Ui(r) 



8-KfJ, 



r j r k Sjk 



r * + ' ^3 



(8) 



with the stresslet tensor S given by 



= / Sifj(s)ds - -± I s k f k (s)ds - n 
ly Jy Jy 



,y i 



U: 



s)nj(s) + uf(s)rii(s) dS(s). 



(9) 



Note that the definition of t he stresslet obtai ned using the single and double layer potentials is the 



same as the one obtained by iBatchelorl (|l970l ). In the following, we will refer to two different kinds 
of swimmers, pushers and pullers, by analogy to a simple case where the swimmer can be replaced 
by a drag-generating center and a thrust-generating center. In that case, all the components of 
Pij = Jy SifjdS are zero except p\\. For a pusher, the thrust generating center (e.g. flagellum) 



is located behind the drag -generating center ( e.g 



configuration and pn > (|Lauga and Powers 



lead) and pn < 0. A puller has the opposite 



20091 ) (note that f was defined as the force density 



from the fluid acting on the swimmer, so a pusher acts with a force distribution on the surrounding 
fluid as directed away from its body along the swimming direction, whereas a puller acts on the 
fluid with a force distribution directed toward the body along the swimming direction). Finally, by 
taking the curl of Eq. (jSJ), it is straigtforward to get that the vorticity field created by the swimmer 



is 



^i(r) 



J a 



4 



(10) 



47r/i r 5 

In this work, we will keep the stresslet tensor S general, to model arbitrary swimming modes. Its 



only constraints are: (1) S T = S in order to enforce torque-free motion, and (2) tr(S) = 0, to ensure 
the conservation of mass through any closed surfaced enclosing the swimmer. In general, S depends 
on the orientation of the swimmer. In the following, we assume that in a frame geometrically 
attached to the swimmer, the stresslet is time-independent in intensity (eigenvalues of the tensor) 
and direction (eigenvectors) so that S = R r 5]R where S is the intrinsic (traceless) stresslet in the 
set of axes B = (e,e',e x e') and R T is the matrix whose columns are the coordinates of B in the 
absolute reference frame Bq. As the swimmer moves, R(t) depends on time but S remains constant. 
Since R is unitary and corresponds to a rotation in three-dimensional space, it corresponds to only 
three degrees of freedom. 

In vector notations, Eqs. (f8j)- (fT0|) become at leading order 



u(r) 



r T • S • r 



8np 

2.3 Coupled motion of two swimmers 

We now consider two identical rotating swimmers, characterized by their position r,- and their 
orientation defined by the two orthogonal vectors ej and (j = 1,2). The corresponding rotation 
matrices Rj are defined as above. The problem is non-dimensionalized using the radius p = Uq/0,q 
of the swimmers' circular trajectory and their intrinsic velocity Uq. The tensors Sj are scaled using 
a particular norm A of Sj (for example the magnitude of its largest eigenvalue) — A is an intrinsic 
property of S and is therefore identical for both swimmers. 

In the far-field approximation, the velocity and rotation of swimmer 1 induced by swimmer 2 are 
respectively equal to the velocity and rotation rate (i.e. half the vorticity) induced by the motion 
of swimmer 2 alone at the posi tion of swimmer 1. We neglect any higher-order term arising from 



the finite size of the swimmers (Kim and Karilla 



199ll ). Such higher order corrections correspond 



to a modification by the presence of swimmer 1 of the velocity field created by swimmer 2. The 
non-dimensional distance r between the two swimmers must therefore satisfy r 3> a/ p. To restrict 
ourselves to the simpler case, we also implicitely assumed that the swimmers are spherical. In the 
case of a non-spherical swimmer, a correction must be added to the rotation rate even in the far- 
field approximation, which physically arises from the alignment of an elongated body in a straining 
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(irrotational) flow (jPedlev and Kesslerl . 119921 ; lLauga and Powersl . 120091 ) . The different limitations 
introduced by these approximations are discussed in §5.21 

Using the results of the previous sections, and under the assumptions presented above, the 
dimensionless equations of motion of the coupled swimmers become 



dri 
~d7 
dr2 

"dT 

dei 
~d7 
de[ 
~dt~ 
de 2 
~dt 
de^ 
dt 



e x - 7 
e 2 - 7 



(ri - r 2 ) • S 2 • (ri - r 2 ) 

|ri - r 2 | 5 
(r 2 - ri) T • Si • (r 2 - n) 



fa - ri), 
x ei, 



x e 1; 



|r 2 - n| 

i , 7(ri ~ r 2) x [S 2 • (n - r 2 )] 

e i + i 15 

|ri - r 2 |° 

7( r i - r 2) x [S 2 • (ri - r 2 )] 
|n - r 2 | 5 

e 2 + 7(r2 -f^^^Ve, 

7(r 2 - n) x [Si • (r 2 - n)] 
|r 2 - ri| 5 



x e 2 , 



(12a) 
(12b) 
(12c) 
(12d) 
(12e) 
(12f) 



where 7 = 3A/(8irij,p 2 Uo). Defining r = r 2 — n, the relative position of the swimmer, and r = |r|, 
their relative distance, these equations can be rewritten for the relative motion of the two coupled 



swimmers as 



dei 
"dT 



de; 



dr 

dt 

7r x (S 2 • r) 
7r x (S 2 • r) 



dt =e2 ~ ei 



x ei, 



7 



r T -(S 2 + Si)-r 



x e 



l) 



de 2 
~dt 
de^ 
dt 



e,+ 



7r x (Si • r) 



x e 2 , 



7r x (Si • r) 



x e 



2' 



(13a) 
(13b) 
(13c) 



Defining ro = (n + r 2 )/2, the position of the midpoint between the two swimmers, the global 
motion of the pair of swimmers is given by 



2— = e 2 + ei 
dt 



7 r T • (Si - S 2 ) • r 



(14) 



Eq. (|13p is a system of five vector equations for r, ej and e'- (j = 1,2), which is closed because 
the knowledge of ej and e'j entirely determines Rj and therefore Sj. Moreover, the equalities 
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• e'j = and ej ■ ej 



e'- ■ e'- 



1 for j = 1,2 mean that, a priori, Eqs. (|13p - (|14p correspond to a 



twelve-dimensional dynamical system. Eq. (|13f) can be solved first for the relative motion since it 
does not involve rg, and one can then obtain the absolute displacement tq from Eq. §H 



3 Far-field interaction of two rotating swimmers 



We are interested in this section in the behavior of Eq. (|13p when the swimmers are far from each 
other, namely, when their relative distance is much greater than the radius of their trajectory 
(r 3> 1). We can focus our attention to the relative motion of the swimmers defined by r as their 
absolute mean displacement ro does not influence Eqs. (|13p . 

Rescaling the distance between the swimmers as r = r*/e with e <C 1 and r* = O(l), the 
equations for the relative motion are obtained from Eq. (113p as (dropping the star superscripts for 
clarity) : 



dr 

dt 
dei 

de 2 



e(e 2 -ei) + e 3 F(r, ei, e[, e 2 , e' 2 ) 



de; 



e[ x ei +e 3 Gi(r,e 2 ,e 2 ) xe b — ^ = e 3 Gi(r,e 2) e' 2 ) x e[, 

de' 

e' 2 x e 2 + e 3 G 2 (r,ei,ei) x e 2 , = e 3 G 2 (r, ei, e' x ) x e' 2 , 



(15a) 
(15b) 
(15c) 



with 



F(r,ei,e / 1 ,e 2 ,e' 2 ) 

Gi(r,e 2 ,e 2 ) 
G 2 (r,ei,e' 1 ) 



r T -(S 1 + S 2 )-r 



7r x (S 2 ■ r) 
7r x (Si • r) 



(16a) 

(16b) 
(16c) 



which are at most O(l). In addition, differentiating the equations for e.j in Eqs. (|15b -c) with respect 
to time leads to 



dt 2 



+ ei 



dt 



x e.; + e J 



dGi 

dt 



x e, + (Gj • e»)e- - 2(e- • G t )e 2 



+ e 6 Gj x (Gj x ej 



(17) 
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since ej • e[ = 0. 



3.1 Multiple-scale analysis 



The equations for e-i in Eq. (|15p suggest that in the limit of small e there are two different time scales: 
The short time-scale is O(l) and corresponds to the intrinsic circular motion of the swimmers, 
whereas the long time scale is 0(e~ 3 ) and corresponds to t he motion induced on one swimmer 



1978 ) with the 



by the other. Using the formalism of multiple-scale analysis (jBender and Orszael . 
assumption of scale separation arising from the far-field approximation (e <C 1), we now formally 
consider all the fields as functions of two variables t and r = e 3 t. The time derivative operator d/dt 
must then be replaced by d/dt + £ 3 <9/<9r, and the different vector fields are obtained as regular 
perturbations series in e 

r = r (°) +er (l) +£ 2 r ( 2 ) +..., (18a) 
e, = e f +ee«+e 2 ef > + ..., (18b) 
e> = ef +e ef) +e 2 ef > + ..., (18c) 

and the functions F, Gj can also be expanded as power series in e, each term being computed from 
the expansion of and r. Introducing this expansion in Eq. f)15|> . we obtain the dynamical system 
at successive orders, which we now solve. 
At order O(l), we have 

at - at - at ~ e ' ' ' 1 s ' 



at order 0(e) 



drW (o) (o) de[ {1) de\ 1} /(0) (1) (0) 

- e> ' x e] ' + e> ' x e\ , (20) 



dt 2 1 ' dt dt 
at order 0(e 2 ) 

5l * (2) (1) U) de 'i 2) n de i 2) '(0) (2) , '(1) (1) , '(2) (0) , 01 s 

-^=e y 2 >-e\>, -^- = 0, -y=e, W xe)'|e 1 u xe; j +e j w xe)', (21) 
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and at order 0(e s ) 



Q r (0) Q r (3) 




(22) 




dr dt 



Note that we are only interested in the leading order behavior of each function. Eq. (I19p gives that 
the leading behavior of r and only varies with the long time scale r. However, it is necessary to 
go up to the terms of order 0(e 3 ) to obtain the r dependence of these functions. This results from 
the ratio between the two time scales being 0(e 3 ) while the first correction to r is 0(e). We note 
from the structure of Eqs. (|20p - (|22p that the i-dependance of the 0(e J ) term in r is determined by 
the previous order in the expansion of e,;. We also note that the relation between and ef ^ is 



linear. 



If we now introduce the expansion from Eq. (|18p into Eq. (|17p . we obtain 



Q2 e (j) 



+ ey> =0, < j < 2. 



(23) 



This equation can be integrated in t as 



ef>{t, t) = a^(r) cos(i) + bf >(t) sin(i). 



(24) 



If we note (.) the t-averaging operator between t and t + 2ir, we therefore obtain 




< j < 2, 



(25) 



and therefore (V 1 ) ^ and (r^ 2 ) ^> are functions of r only. We can now take the average of the first 
equation in Eq. f)22[) remembering that r^ ^ is a function of r only 



dr(°) 



(F(r W , e( Q) , ef) , ef , e? 0) ) ) = r ( 3 > (t , r ) - r ^ (t + 2vr, r) 



(26) 



dr 



13 



Prom Eqs. (|19p - (|2ip . is independent of t for < j < 2. Therefore we have 

x e[ {j) = ap } (r) cos(i) + } (r) sin(t). (27) 

From the definition of F and Eqs. (|24p - (|27p . we can write 

F(r(°\e¥\e} \dp,f4 0) ) = A(r) cos(2t) + B(r) sin(2t) (28) 

+C(r) cos(i) + D(r) sin(i) + E(r), 

and the left-hand side of Eq. (|26p is a function cx(t) of r only. Then we have r^ 3 -*(t + 2ra7r,r) = 
r( 3 )(t, r) — na(r). For the perturbation expansion assumption to remain valid, a must be equal to 
zero. Therefore, 




The same procedure applied to the second equation in Eq. ([22]) gives 

ef>(t + 2vr,r) - ef (t,r) = + (g^, ef , ef )) . (30) 

From the definition of G%, Gj(r^ \ ) can be written in a similar form as F(r^, , , , ) 
in Eq. (|28p . The right hand side of Eq. (|30p is therefore a function of r only and to avoid secular 
terms, both sides of the equation must be zero and 




At leading order, the system behaves therefore as 

r = r(°)(r) + 0( e ), e t = e^^r) + 0(e), ej = ef°\r) + 0(e), (32) 
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Figure 2: Multiple-scale analysis for the motion of the two swimming cells: The leading order 
motion is characterized by the distance between the mean positions of the two swimmers r and the 
orientation of their rotation vectors e' x and e' 2 . These three vectors evolve with the slow time scale 
r, while the instantaneous position of each swimmer is the superposition of their mean and relative 
motion on the slow time scale r and the circular motion on the fast time scale t. 



with 




(33a) 
(33b) 



The different notations are summarized on Fig. [21 where the superscript (0) was dropped for 
clarity. We note that to achieve our final result, the hypothesis ej • = was crucial: It is only 
because the intrinsic motion produces no net displacement over a period that the separation of 
scales is possible. If it is not the case but the dot product of these vectors is small, the intrinsic 
trajectory would be an helix but the net displacement h over one period would still be small. 
We expect that the analysis remain valid provided h <C r, but this should be confirmed with a 
perturbation expansion in the helix step, which gives a new small parameter. 

3.2 Computation of the average quantities 

In this section, we compute quantities such as (f e'j°\ e^, e'^ and (^Gi(A°\e^ , e^ "*)^) 
with the average taken over one period of the short time-scale t. For clarity of notations, we drop 
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Figure 3: Notations for the computation of the average quantities (F(r, ei, e^, &2, e^)) and 
/Gj(r, ej, over a period of the short time scale t corresponding to one period of the circu- 
lar motion of swimmer 1. The vectors in black are constant over this time-scale (they depend on 
r) and the grey vector ei evolves as Eq. (|34|) . 

the (0) exponents with the understanding that we are only considering vector fields of that or- 
der. Over this period, r and e^- are constant vectors. Defining a unit vector i orthogonal to e[ 
and r, the basis B p = (i, x i,eQ is orthonormal (Fig. [3|). The instantaneous intrinsic directions 
corresponding to the intrinsic translation and rotation velocities vary as 

ej = cos ti + sin t e' { x i, x e, = - sin ti + cos te[ x i (34) 

with no loss of generality since we can redefine the origin of time so that ej is orthogonal to r at 
t = (r • i = 0). The vector r can also be decomposed in B p 

r = rie^ + ^e^ x i with r\ + r\ = r 2 . (35) 
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Then we easily obtain 



R r 



e, • r 



e, • r 



( T2 sin t 



ri 

-J"2 cos t 



and 



<(e 4 -r) 2 > = ^, ((e^-r) 2 > = r 1 2 , (((ejxe^rf)^, 
((e 4 • r)(e< • r)> = ((e, • r)[(e' t x e,) • r]> = ([(ej x e,) • r](ej • r)> = 



and therefore 



(r T -S 4 -r> = ((R-r) T -S-(R-r)> 



(£ n + £ 33 ) + r 2 £ 
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[r 2 - (e; • r) 2 ]tr(£) 



Finally, since tr(S) = 0, 



(o) jo) '(0) p (0) /(0)A _ 7S 



F(r^,er,er,er,er) 



Similarly, 



and 



-22 



2r 2 - 3 [(ei • r) 2 + (e' 2 • r) 



Si • r = [Snr 2 sin t + Si 2 ri - Si 3 r 2 cost] e, 
+ [E21T2 shit + S 2 2n - S 2 3?"2 cost] e- 
+ [E 3 ir 2 sint + £ 32 ri - £ 33 r 2 cost] (e* x e-J 



e, x r = n sin t i — n cos t (e£ x i) + r2 cos t e£, 



e„ x r 



-r 2 i, 



(e; x ej) x r = — n cos t i — r\ sin t (e^ x i) + r2 sin t e[ 
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from which we obtain after time-averaging, 



<(S<T)xr> 



rir 2 



(En — 2£ 22 + S33) i + (£31 — £13) 



2„/ 



nr 2 , . 
e- x 1 



(43) 



The last term in the last equation is equal to zero as £ is symmetric, and identifying r\ = • r 
and r 2 i = r x e^, the previous equation becomes: 



((Si-r)xr) 



35^22 / / 



e'i-r)(r x ej). 



(44) 



Therefore, 



Gi(r 



(0) (o) /(o) 

) e 2 > e 2 



G 2 (r(°),eS ),ef ) 



3 7 E 22 (e 2 -r)(r xe 2 ) 
2r 5 

37S 2 2(e / 1 T)(rxej) 
2r 5 



(45) 



Finally, the relative equations of motion for the slow varying fields r, and e 2 become with 

A* = 7 S 22 



dr 
dr" 



3- = M 



2r 2 - 3 [(ei • r) 2 + (e' 2 • if 



de^ 
~d7 
de^ 
dr 



2r 5 

3^(e 2 -r)[(r x e 2 ) x e[] 

2r 5 

2r 5 



(46a) 

(46b) 
(46c) 



The absolute motion of the swimmers can be determined on the long time scale r by averaging (|14p 
over the short-time scale and obtain 



d(r ) 3// 



dr 



^[(e 2 . r ) 2 -(e;.r) 2 ]^ 



(47) 



A comparison of the dynamical systems given by Eq. (|46|) and Eq. (|13|) shows that the averaged 
equations, Eq. (06]), correspond to the interaction of two stresslets of equal intensity 3[/,/2(e' 1 ei—l/3) 
and 3^/2(e 2 e 2 — 1/3) respectively located at the mean position of swimmers 1 and 2 with no intrinsic 
velocity. This suggests that a single swimmer creates an average far-field in the form of a stresslet 
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whose intensity is 3/i/2 and whose orientation is entirely determined by its intrinsic rotation vector 
e^. This statement is proven rigorously in §3.31 The intensity of the averaged stresslet is equal to 
3X22/2, where X22 is the diagonal component of the instantaneous stresslet along the direction e^. 
We observe that all the other components of S disappear in the averaging process. 

By analogy with the case where the instantaneous stresslet is equal to a force dipole, resulting 
from the superposition of a drag force and a thrust force, we will consider in the following two 
kinds of swimmers: 

• Pushers with fi > 0: In this case the thrust generating center is located behind the drag 
generating center; 7X11 < and 7X22 = 7X33 > with all other components equal to zero 
[see Eq. flUJ)]. This is for example the case of a swimmer with a nagellum located behind its 
drag-generating head, such as spermatozoa, or most flagellated bacteria. 

• Pullers with fj, < 0: In that case, the thrust is generated in front of the drag-generating center; 
7X11 > and 7X22 = 7X33 < [see Eq. Q]. This is for example the case for swimmers using 
their flagella in a breaststroke pattern to pull their bodies, such as the alga Chlamydomonas. 

It is important to point out here that we manage to obtain a system of equations for e' l5 e 2 and r 
only, but that the position of each swimmer on its instantaneous circular trajectory is not important 
— in particular the relative phase of these instantaneous motions. Two conditions are necessary for 
this simplification to occur. First, the average flow field created by an isolated rotating swimmer is 
independent of time and also independent of the direction of motion on the circular trajectory (see 
the following section). This is a consequence of the fact that the instantaneous flow field created by 
the swimmer does not have any azimuthal component. The second condition is that the swimmers 
are spherical, and the averaged velocity induced on swimmer 2 by swimmer 1 only depends on the 
properties of the averaged flow field induced by swimmer 1 and not the orientation of swimmer 
2. This would not be the case if the swimmers were non-spherical: then, the induced velocity and 
rotation created by swimmer 1 on swimmer 2 would not only depend on the position and trajectory 
of swimmer 1, but also on the orientation of swimmer 2 with respect to the principal axes of strain 
of the local flow (see the discussion in §5.2|) . For non-spherical swimmers, the averaging process is 
more subtle and the phase of the instantaneous motions of the two swimmers does not disappear 
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in the averaged equations; it remains however a constant parameter of the problem since both 
swimmers have the same intrinsic translation and rotation velocities. 



3.3 Far-field averaged velocity field created by a rotating swimmer 

The results of the previous section suggest that, on average, a rotating swimmer behaves like a 
stresslet in the far-field. We explore this result in more detail in this section. The behavior of 
the far-fiel d velocity i s of in terest to characterize the rheological properties of a suspension of such 



swimmers (Batchelor 



19701 ) . In this section only, we consider an isolated swimmer, and compute 
the time-averaged flow in the far field. The swimmer trajectory is a circle oriented by its rotation 
vector e' parallel to the vertical axis and we choose the origin of the reference axes as the average 
position of this swimmer. Let denote by e(t) the instantaneous position of the swimmer (|e(t)| = 1 
by our choice of scaling) and e its velocity vector. If i is an arbitrary constant unit vector orthogonal 
to e', we can define the origin of time such that: 



e(t) = cos ti + sin t e' x i, e = — sin t i + cos t e' X i. 



(48) 



We are interested in the velocity field created by this swimmer at a position x far from the origin 
(i > 1). The instantaneous velocity field at x is given from Eq. (fTTj) by 



u(x) 



"7 



(R • r) T • £ • (R • r) 



r, with r 



(49) 



and 



R r 



e • r 



e ■ r 



e • x 



e • x 



e • e 



e ■ e 



e • x 



e • x 



e x e 



■r J 



\ (exe')'X y 



\ (e x e') • e y 



\ (exe')-x j 



w 

(50) 
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Therefore from Eq. (|48|) . noting once again (.) the averaging operator over a 2-7r-period, we have 



<^ 2 > 

{Pi) 
(Pi) 
(P2P3} 
(P1P2) 



l r 
2 



(x-(e'xi)) 2 + (i-xr 
,2 



(x.e') ; 



1 

1+ 2L 
-xe' 



(P1P3) = 0. 



x 2 - fx-eV 



(51a) 
(51b) 
(51c) 
(51d) 
(51e) 



Keeping only the dominant terms, we have on average 



<r T • S • r) 



S11 + S 



(x.e') : 



+ S 22 (x.ef = ^[3(x.e') 2 -x 2 



We also have 



1 1 



e • x 



— = — 1 + n^- + - 



(52) 



(53) 



Since we are interested only in the dominant term in the far-field averaged behavior, we write 



(R • r) T • £ • (R • r) 



(Rt) t -E.(Rt) 



X" 



(54) 



as all the corrections to this expression are of higher order in 1/x. Grouping all terms, we finally 
obtain the far-field averaged flow 



(u) (x) 



7S22 [x r • (3eV - I) • x] x 



(55) 



We recognize here the velocity field created by a steady stresslet 3///2 (e'e' — 1/3) consistently with 
the results of the previous section. Physically, the results of Eq. (|55p indicate that, for cells which 
behave instantaneously as pushers (pullers), the averaged flow is that of a puller (pusher) along the 
axis of rotation of the circular motion. 

We observe in Eq. (|55p that the average flow remains identical by changing e' into — e': the 
average flow is therefore not modified by a reversal of the circular motion (along the same trajec- 
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tory). 



4 Analysis of the far-field interaction 

4.1 Reduced forms of the equations 

We now return to the coupled equations derived using the multiple-scale analysis. Defining the unit 
vector e z of the direction between swimmer 1 and swimmer 2, e z = r/|r| and by differentiation in 
time we obtain 

dp. 1 rl r I v rlr\ 

(56) 



dr 



r. 



\r\ dr \ | r | 3 dr J 

But from Eq. (|46h . we note that dr/dr = lZr, with 1Z a scalar function of r and e^. Using this result 
in (|56p . we obtain that e 2 = r/|r| is a time-independent unit vector set by the initial conditions. 
The mean distance between the two swimmers maintain a fixed direction. In the following, & z 
denotes the fixed direction between the two swimmers' positions. The vectors e[ are defined from 
e z by their polar and azimuthal angle 0, and fa. Choosing two constant unit vectors e x and e y so 
that ( e x ,ey,e z ) is orthonormal, then 



= sin 6i cos fae x + sin Oi sin fae v + cos 9i& z . 



(57) 



Note here, that the definition of fa depends on the definition of and e y which can be rotated 
arbitrarily in the plane orthogonal to e z . Therefore, only the intrinsic £ = fa — fa has a physical 
meaning. Then in the frame (e x ,e y ,e z ) we have 



and 



/ 



(ei • r)[r x e[] x e' 2 



cos Wi cos f 2 sin f i cos <p\ 
cos d\ cos 02 sin 6\ sin fa 
cos#i sin^i sin 6*2 cos{fa — fa) j 



de 2 
dr" 



d02 
dr 



COS 02 COS fa 



cos 02 sm <p2 
— sin 62 



+ 



dr 



' — sin 02 sin fa 



sm y 2 cos <p2 




(58) 



(59) 
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By identification, the system given by Eq. (|46p can then be rewritten as a four-dimensional dynam- 
ical system 



^ = A [2 - 3(cos 2 6 X + cos 2 2 )1 , (60a) 
dr zr z L J 

— - = — ^ cos #2 sin 6*2 cos £, (60b) 
dr 2r d 

d6l 2 3 ^ a • a t fan \ 

— — = — ? cos #i sin Ui cos £ , (60c) 

dr 2r d 

sin sin ^2 -r~ = ^ cos 6q cos #2(sin 2 6*i + sin 2 #2) sin£, (60d) 

dr 2r 6 

where we have used £ = fa — 4>i- The notations for Eq. (I60j) are summarized on Fig. [H Note that 
Eq. (f60l) can be simplified even further by defining a = 2r 3 /3/U, Xi = cos #j and y = sin#i sin 02 cos£, 
and we obtain 



da 
dr" 



3(xf + sl), (61a) 



dxi 

a—— = -x 2 y, (61b) 
dr 

drc2 

a—— = -xiy, (61c) 
dr 

a— = x\X2{2 — x\ — x\) . (61d) 
dr 

Physically, a is proportional to the third power of the distance between the swimmers. It is negative 
for fj, < (pullers) and positive for \x > (pushers). From the original physical problem, we also 
have the following three mathematical constraints: 

• The variable a is either positive or negative. A change of sign of a requires a cancellation 
of r at a finite time and a collision of the swimmers. Such a collision obviously violates 
the far-field approximation, and the present theory is not valid when a gets small. In the 
following, we will refer as "collisions" to regimes where the present theory predicts a decrease 
of the relative distance to an arbitrary small number, at which point additional modeling is 
required. We will therefore focus on solutions for which the sign of a is fixed. 

• The variables x\ and X2 are cosines, therefore —1 < {x\,X2} < 1. 

• From the definition of y, < y 2 < (1 — x\){l — x 2 ). 
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Figure 4: Definitions of the various variables for the average motion (see text for details). 
4.2 Relative equilibria and stability 

We focus here on relative equilibrium positions, for which on the long time scale, the swimmers 
do not move relatively to each other. There can be however a mutual motion of the swimmers 
(d(r )/dr 0). 

4.2.1 Equilibrium points 

From (|6ip . there is only one type of equilibrium points obtained for (a, x\,X2,y) = (ao> ±a/2/3, 0, 0), 
or symmetrically (a, x\, X2, y) = (acb 0, ±\/2/3, 0) for any value «o of a. 

Physically, the distance between the swimmers can take an arbitrary value but the orientations 
of the rotation vectors must correspond to a very specific configuration. One swimmer's rotation 
axis makes an angle cos" 1 a/2/3 ~ 35 ° with the distance between the swimmers. The second 
swimmer's rotation axis is orthogonal to the first swimmer's and their relative distance (so e^- is 
orthogonal to the plane defined by and e z , with j ^ %). 
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The linearized system about one such equilibrium is obtained as: 
/ ™ _ ™„ \ / n —.Hk n n \ 



d_ 

dt 



V 



a — «o ^ 


/ 





-V6 


x i - \/i 




















f / 


V 















O 



4 

3ao 







' a — ao 



/ 



V 



x 2 
V 



(62) 



The eigenvalues of the above matrix are A = ±i2\/2/3ao, and A = with multiplicity 2. The 
dimension of the subspace associated with A = is however equal to 1 . It is therefore not possible 
to conclude from the linearized system on the stability of the equilibrium of the non -linear system 



as one of the eigenvalue of the linearized system is identically zero (neutral stability) ((Sastryj, 
We will show rigorously in §4.41 that this equilibrium is unstable. 



19991 ) 



4.2.2 Rotational equilibria 

Another situation of interest is the case where the direction of the circular motions, remains 
fixed relatively to e 2 . Only a (or equivalently the distance between the two swimmers) depends on 
time. This occurs for two different configurations. 



Swimmers with same axis of rotation: x\ = ±1 and x 2 = ±1. The two swimmers have quasi- 
circular trajectories in two parallel planes and e[ are both aligned with e z . As a direct consequence 
of the definitions of Xi and y, y must be zero at all time. From Eq. (|61ap . the evolution of a in 
that configuration can be computed 



a = a = «o — 4t, and r = r 



(63) 



The overbar denotes the reference configuration (rotational equilibrium) we are considering. There- 
fore, x\ = ±1 and X2 = ±1. Swimmers with [i < (pullers) tend to repel each other while swimmers 
with [i > (pushers) attract each other, until the scale-separation assumptions of the multiple- 
scale analysis break down. We observe that the collision time scales like r\j \i ~ Ft? / pal with R the 
dimensional distance between the swimmers, p the radius of their circular trajectory and / and a 
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the length and head radius of the swimmer, respectively. 

The stability of this time- varying configuration is now investigated by decomposing each variable 
/ (with / = a, x\, X2, y) as / = / + /' and /' is a small perturbation. At leading order, Eq. ([ST]) 
can be rewritten 



a 



dx\_ 



da' . _ . 

— = -Q(x 1 x 1 + x 2 x 2 j 



-X2V 

dV 

dr 



- dXn 

a— — = -xiy 
dr 

— 2X1X2 — X2X 2 ) 



(64a) 
(64b) 
(64c) 



or equivalently 



da' 

— — = -6(xix 1 + x 2 x 2 
ar 



2J> 



_ d 
a— 
dt 



( J \ 



\ V J 



V 









-X2 








-XI 


2x 2 


-2xi 






7 V y ! 



(65) 



where the bar quantities correspond to the relative equilibrium (a = uq — 4r and x\ = 1) and the 
prime quantities are perturbations. 

The last system can be solved exactly if diagonalized. Defining 



/ 



•?'2 
X\ 





Z2 

\ Z3 ) V 

it decomposes into 



-Xl Xl 

-x 2 x 2 
2 2 



\ 



x 



x 2 



V y / 



2x 2 

Xl 
-Xl 



-2xi 
x 2 1 
—Xi 1 



\ 



( Z X \ 



, (66) 



dgi _ d£2 
dr dT 



2z 2 dz 3 2z 3 



ao — 4r dr ao — 4r 



(67) 



which can be integrated easily into 



Z\ = Zl,o, Z 2 = Z2,0 1 , ^3 = ^3,0 1 

a J \ a J 



(68) 



The original variables are obtained by linear combinations of these solutions, and we observe that 
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the configuration is unstable with algebraically growing perturbations. 

Figure [5] illustrates this situation and compares the prediction of the far-field model for the 
averaged motion (Eq. [B"T|) to the full set of equations (Eq. [T3j) . Considering two pushers (/j, > 0) 
that have initially almost the same axis of rotation {61,62 <C 1), the hydrodynamic interactions 
create a mutual attraction of the swimmers following the approximate law (Eq. [63]). This rotational 
equilibrium is unstable and as they get closer from each other, the planes of the trajectories of the 
two swimmers undergo a quick rotation, bringing the two swimmers from a co-axial to a co-planar 
configuration (see next section) in which the interaction of the two pushers have now a repulsive 
effect. Figure [5] also allows to show the agreement between the simplified model (Eq. [6T]) and the 
full equations of the system. 



Two-dimensional configuration: x\ = x 2 = and y = ±1. For 6i = ir/2 and £ = (y = 1) 

for co-rotating and £ = tt {y = —1) for counter-rotating swimmers, both swimmers are in the same 
plane with their rotation axes orthogonal to the plane of motion. Note that the two-dimensional 
configurations are actually only particular cases of orientational equilibria: X\ = and — 1 < y < 1 . 
As above, this configuration is a rotational equilibrium only, as the distance between the swimmers 
varies in time according to 



1 /3 

a = a + 2r, r = (rjj + 3/ir) 



(69) 



This time, swimmers with /i < (pullers) attract each other while swimmers with fi > (pushers) 
repel each other. Here again, a stability analysis can be performed, and the linearized dynamics 
becomes 



d Q ' n( '2 1 '2\ 

— j— — 0\X\ T X 2 )i 



which can be integrated exactly as 



_ d 
a— 
dt 
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(70) 



2/o > 



30i ~\~ Xo 



1 + — 

a J 



1 + — 

a J 



(71) 
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Figure 5: Interaction of two pushers (/i = 1) with circular motions that are initially almost coaxial 
(di = 02 = 0.22, £ = 0.68) and an initial distance equal to r = 11.4. (a) Trajectories of the two 
swimmers: the initially coaxial pushers attract each other until hydrodynamic interactions modify 
the orientation of their circular trajectories and they become co-planar, leading to a repulsive 
interaction, (b) Evolution of the distance between the swimmers and (c) evolution of the parameter 
a. In both (b) and (c), the light grey curve corresponds to the full equations (Eq. [13]) for which 
the circular motion of each swimmer is resolved, and the black curve corresponds to the simplified 
model for the averaged motion (Eq. IBTI) . Note that the two curves agree with each other until the 
swimmers get close to each other. 
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once again leading to the instability of these configurations with algebraically growing perturba- 
tions. 

4.3 Reduction to a two-dimensional problem 
4.3.1 Conserved quantities 

The system given by Eq. (|6ip can be simplified even further by observing that 

A = x\ - x\ = cos 2 6»i - cos 2 9 2 (72) 

is a conserved quantity. Without any loss of generality, we can assume A to be positive (the 
equations are symmetric with respect to a switch between x\ and x^). In the (x\, X2)-plane, the 
system moves along a hyperbola, and we can define the parametric coordinate a such that 

x\ = \J~A~ cosher, x 2 = v^Asinhu. (73) 

To be rigorous, x\ should be equal to ±t/A cosher. However, one can change (x\,X2,y) into 
(— x\,— X2,y) by changing the definition of e z to — e z (or equivalently, switching the indices of 
the swimmers), and we therefore restrict ourselves to xi > by redefining er appropriately. Intro- 
ducing this change of variables into Eq. (|61cj) leads to 

y = —ad. (74) 

Using this relation in Eq. (|61d|) . we obtain 

— a(aa + aa) = — sinh 2<r (2 — A cosh 2a) , (75) 

and multiplying by & and integrating with respect to time, we obtain 

A 2 

a 2 a 2 = -Acosh2cr + — cosh 2 2cr + C, (76) 
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where C is a constant of integration. Noting from Eq. (|73p that A cosh 2a = x\ + x|, we have 
therefore proven that 



A = *?-*i and C = y 2 + (* 1 + ^)(l-^I±^) 



are two conserved quantities in this problem. Finally, defining the new variable X = x\ + x\ 



j4cosh2cr, the system given by Eq. (|61[) is equivalent to 



aX = e^{X 2 - A 2 ) (4C - 4X + X 2 ) (78a) 
d = 2-3X (78b) 

with e = ±1. We have therefore transformed the four-dimensional system, Eq. (|61|) . into a two- 
dimensional system, Eq. (|78p . The values of the constants A and C, as well as the initial values of 
X and a can be obtained from the initial conditions of the four variables (a, Q\, 02, £)• The choice 
of the sign of e is discussed in §4.3.31 

4.3.2 Bounds on the different variables 

From the constraints detailed at the end of ^4.11 and the definitions of the variable X, and the 
constants A and C, we have the following four constraints. 

• < {xfjxl} < 1 therefore 

< A < 1, (79a) 
A < X < 2 - A. (79b) 

• From Eq. (|77|> . C = v 2 + X-X 2 /A. Using the previous bounds on X as well as the inequality 
y 2 < (1 — x\){\ — x 2 ), we obtain that 

C<l-~ (80) 
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Prom Eq. ([77| . we have y 2 = C - X + X 2 /4, and therefore 

X < 2 (l - Vl^C^j . (81) 

Because of Eq. (|5U|) , we have 



2-A>2 (l-VT^C , (82) 



and Eq. (|8ip is actually a tigher upper bound than Eq. (|79b ). 

• Finally, Eq. ([ST]) and X > A implies that C > A - A 2 /A. 

In summary, the following inequalities must be satisfied 

A 2 A 2 / \ 

0<A<1, A- — <C <1-—, A< X < 2(1- VT^C) . (83) 

4.3.3 Choosing the sign of e 

For given values of A and C, and given initial conditions on X and a, there are two possible 
solutions corresponding to e = ±1 initially. In Eq. (JTBaJ), the square-root of the right-hand-side is 
positive and e has therefore the sign of aX = —AxiX2y. Differentiating with respect to time and 
using Eq. (|6T|) . we obtain 

and X and X are continuous functions of time. 

From the constraints of §4.11 we are only interested in solutions where the sign of a is fixed. 
Therefore, we are not interested in the solutions of Eq. (|78p past a zero of a. The left-hand side of 
Eq. (JTBaj) vanishes only for vanishing X or for collisions. We prove here, that if at t = to, aX = 0, 
then e must change of sign at t = to if a(to) ^ 0. Such a cancellation of the left-hand sign of 
Eq. (|78ap happens only in two configurations 

1. X = A = X m i n or equivalently X2 = (one swimmer's rotation axis is orthogonal to the 
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distance between the two swimmers). For X to reach a minimum at t = to, X < for t < to 
and 6(£q) = e~ = — sgn(a(to)). For small |i — to|, we obtain using Taylor expansion and 
Eq. flgj" 



aX 



4A /„ A 



. C- A+ ^ 
a V 4 



(t-t ), (85) 



which is positive for t > to, therefore e(i^) = e + = sgn(a(to)) = — e ■ 

2. X = 2(1 — \/l — C) = ^mai or equivalently y = 0. For t < to, we therefore have e" 
sgn(a(to))- For small |t — to|, we obtain 



2 v / l _r C (4(1 - \/r^C) 2 - A 2 ) 



a 



(t-to), (86) 



and e + = — sgn(a(to)) = — e . 

With the analysis above, we see that for given values of A and C the system can be represented 
solely in the (X, a) plane. However, if one wants to look at maps of the flow, two maps should be 
superimposed e = 1 and e = —1, one for trajectories of decreasing X and the other for trajectories 
of increasing X. 

4.4 Possible regimes in the far-field interaction of two rotating swimmers 
4.4.1 Monotonic variations of a 

From Eq. (|78b[) . we see that a is an increasing (decreasing) function of time if X < 2/3 (X > 2/3). 
If 2/3 is out of the bounds imposed on X by Eq. (|83p . a and the distance between the swimmers 
are monotonic functions of time. Two such cases can occur. 

1. If A > 2/3, then a < 2 - 3A < 0, and 

• if ao < 0, a ^ — oo and the swimmers get further and further away from each other, 

• if instead a$ > 0, a — > and a collision occur at a finite time (since the time derivative 
of a is negative and has a non-zero negative upper bound. 



2. If C < 5/9, then a > 2 - 6(1 - VT^) > and 
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• if ao < 0, a — > and a collision occurs at finite time (since the time derivative of a is 
positive and has a non-zero positive lower bound), 

• if instead ao > 0, a — > oo and the distance between the swimmers is unbounded. 
4.4.2 General case: Theory 

If A < 2/3 and C > 5/9, then we can prove that X oscillates from its lower bound X m i n = A to its 
upper bound X max = 2(1 — y/1 — C). This statement could be proven rigorously from the equations 
for X and a. We only provide here a qualitative argument for clarity. Since X only vanishes at 
these bounds, X varies monotonically from one to the other. If X doesn't reach the next bound 
(even as t — > oo), then it would have a finite limit and X must go to zero as t — > oo while aX 
remains finite; this combination can only occur if a is unbounded. Therefore, X oscillates between 
its bounds unless \a\ — > oo. 



Then, let t n be the successive times at which X reaches either X m i n or X max and the correspond- 
ing values a n . We are interested in the gain G n = \a n+ \/a n \ and the time interval r n = t n+ \ — t n 
between two sign reversals of X. Prom Eq. (I78p . we obtain that, over an interval where X has a 
given sign, we have 



x 



(2 - 3X)dX 



F(X;A,C,Xo) = y ' = = elog 

1 ^ Jxo V / (X 2 -A 2 )(4C-4X + X 2 ) 



a 



on 



(87) 



T is well defined for X m i n < X < X max as the singularities at the end points are integrable. Using 
Eq. ([8"7]) between t n and i n +i> we obtain the following. 

• If ao > (therefore a > at all time at least until collision), then when X varies from X m i n 
to Xmax, e is positive, and 



logG n = g(A,C), with g(A,C)= / 



(2 - 3X)dX 



x min v / (X 2 -A2)(4C7-4X + X2) 
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• If ao < (therefore a < 0) , then when X varies from X m i n to X max , e is negative and 

logG n = -g(A,C). (89) 

We note here that G n is a function of ^4 and C only and therefore not a function of n or a n . We can 
therefore summarize these results for all initial choice of (a, xi,X2,y) or equivalently (ao, Xq, A, C): 

1. If ao < (puller) and Q(A, C) > 0, a — > and there is a collision between the swimmers. 

2. If a < (puller) and G(A, C) < 0, a — ► — oo and the swimmers move away from each other. 

3. If ao > (pusher) and Q(A, C) > 0, a — ► oo and the swimmers move away from each other. 

4. If ao > (pusher) and G(A, C) < 0, a — ► and there is a collision between the swimmers. 

Note that the particular cases discussed in the previous section (A > 2/3 and C < 5/9) are also 
included in this analysis: the integrand in Q has then a fixed sign. 

It is important to point out that Q determines the regime (divergence or collision) of the two 
swimmers and is a function of A and C only The regime is therefore entirely determined by these 
two quantities, which are only functions of the relative orientation of the rotation vectors of the 
swimmers and independent of their initial separation distance. The maps of the regimes obtained 
for pushers (fi > 0) and pullers (/x < 0) are displayed on Fig. [6J Note that at the boundary between 
the collision and divergence domains, we have Q = 0: The distance between the swimmers remains 
unchanged between t n and £ n +2 and the motion is periodic. This corresponds to a limit cycle. The 
boundary between the regimes shown on Fig. [6] (solid line between the divergence and collision 
regions) can be obtained numerically by finding the values of A and C for which Q(A, C) = 0. The 
other two solid boundaries correspond to C = 1 — A 2 /4 and C = A — A 2 /4 [see the constraints on 
C in (Ell. 



4.4.3 General case: Numerical simulations 

As a followup to our theoretical analysis, we illustrate here the four different possible regimes 
obtained in §4.4.21 These results are displayed in Fig. [7l If A > 2/3 or C < 5/9, the variation of a 
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(a) Pushers (jj, > 0) (b) Pullers (p < 0) 

Figure 6: Maps of the general regime in the (A, C)-plane for the long-time relative behavior of two 
swimmers with positive fx (left, pushers) or negative \i (right, pullers). Note that one can be deduced 
from the other by symmetry, i.e. by changing collision (divergence) by divergence (collision). On 
the left map (pushers), the position of the four examples (a)-(d) of Fig. [7] are indicated. 



with time is monotonic, and can either be divergent (Fig. [Th.) or convergent (Fig. [7b 1 ). If A < 2/3 
and C > 5/9, then the variation of a over half an oscillation is not monotonic: a changes of sign 
when X = 2/3 corresponding to a minimum or maximum distance between the swimmers. This 
leads to a spiral shape of the trajectory in the plane (X, a), and non-monotonic divergence (Fig. [7b) 
or convergence (Fig. [TH) of the relative position between the swimmers. 

4.4.4 Finite time of collision 

We show here that in cases (1) and (4) discussed at the end of §4.4.21 the collision between the two 
swimmers occurs in a finite time. We consider case (4) for example. The time interval between two 
zeros of X is given by 

f Xmax a(X)dX 

T n = t n+i -t n = y ' < a n T(A, C), (90) 

Jx min ^(X*-A>){AC-AX + X>) ~ n K h ^ } 

with 



x mm V(X 2 - A 2 )(AC - 4X + X 2 ) 
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(a) A = 0.1, C = 0.3, X = 0.25 and a = 1 
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(d) A = 0.6, C = 0.7, X = 0.75 and a = 10 

Figure 7: Illustration of the four possible trajectories in the phase plane (X, a) (left), evolution 
of X (center) and evolution of a for four different configurations: (a) monotonic divergence, (b) 
monotonic convergence, (c) non-monotonic divergence and (d) non-monotonic convergence of the 
swimmers. 



3G 



Then we have 



n 



1 _ e (r»+l)6 



T 



1-eO' 



(92) 



fc=0 



as n -> oo and the collision a n = happens at a finite time since £/ is negative. A bound on the 
finite collision time can be obtained in the same way for case (1). 

4.4.5 Analysis of the system equilibrium 

Finally, we know from §4.21 that the system has only one type of possible equilibrium. Using the 
notation from the current section, it corresponds to A = X = 2/3 and C = 5/9, a being arbitrary. 
This point is on the boundary between the collision and divergence domains on Fig. as well as 
on the boundary C = A — A 2 /4. One sees easily from Fig. [6j that for all values of a and for any 
value of the perturbation that does not leave A and C both unchanged, the system will move away 
from its equilibrium. This equilibrium is therefore nonlinearly unstable. 

4.4.6 Absolute displacement of the swimmers through hydrodynamic interactions 

In the previous sections, we have focused mostly on the relative motion of the two swimmers. The 
absolute motion is characterized by the evolution of ro, defined as the middle point between the two 
swimmers. In the far-field approximation, we have computed the average velocity of this middle 
point as vq in (f47l) . Using the notations defined in sections 14.11 and [472], vq becomes 



where e z is a constant unit vector giving the direction of the relative distance between the swimmers. 
A = cos 2 02 — cos 2 9i was shown to be a conserved quantity. The absolute motion of the swimmers 
therefore occurs along the same direction as the relative motion, and vo does not change sign. In 
section 14.3-H we relabeled the swimmers so that A is a positive quantity. With this relabeling, the 
net motion of the swimmers occurs along e z in the direction of swimmer 1 for pushers (// > 0) and 
in the direction of swimmer 2 for pullers (// < 0). The net displacement velocity | vol scales like 
1/r 2 , as expected from dipolar hydrodynamic interactions. 



v = - 



3^ A 

P 

4 r 2 
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5 Discussion 



5.1 Summary of results and biological relevance 

The work in this paper focuses on the hydrodynamic interaction of N = 2 swimmers with circular 



traje ctories, such as the spermatozoa of some marine invertebrates ([Goldstein! . 



1977 



Riedel et al 



20051) . This particular situation represents the simplest configuration in which the effects of hy- 
drodynamic interactions can be studied without considering the full iV-body problem with a large 
number of organisms. Indeed, the confinement of the individual trajectories allows the two swim- 
ming organisms to interact on a much longer time-scale than if they were swimming along straight 
lines. The two cells are assumed here to be spherical and identical, but the description of their swim- 
ming stroke is general, and the far-field interaction analysis is valid for an arbitrary stresslet tensor 
(i.e. an arbitrary force distribution at the swimmer surface). In general, the relative dynamics of 
the two cells is described by a dynamical system with nine degrees of freedom. 

In the far-field assumption, a separation of time scales occurs between the period of the intrinsic 
circular motion of the swimmers and the time over which hydrodynamic interactions have an order- 
one effect on their trajectories. As a result, the dynamical system is investigated using a multiple- 
scale analysis. In particular, the average motion resulting from the instantaneous interaction of 
the two swimmers is found to be strictly equivalent to the interaction of two modified stresslets, 
obtained as the stresslet for each swimmer averaged over a period of its intrinsic motion (in other 
words, the time-averaged interaction between the swimmers is equal to the interactions between 
the time-averaged swimmers). Furthermore, the direction of the relative distance between the two 
swimmers is found to be independent of time, and the average problem was reduced to a four- 
dimensional dynamical system for the distance between the swimmers and the relative orientations 
of their rotation vectors. 

We then proceed to a detailed mathematical analysis of the dynamical system. We show the 
existence of one type of equilibrium, which is linearly neutrally stable but nonlinearly unstable, and 
two types of rotational equilibria, which are both linearly unstable with algebraic growth. We then 
show the existence of two conserved quantities, thereby allowing a reduction to a two-dimensional 
dynamical system. We proceed to identify geometrical bounds on the dynamics, and we show that 
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only two general long-time behaviors are possible: Either the swimming cells swim away from each 
other, or they get closer from each other (until the far- field assumption breaks down). In these 
divergence and collision scenarios, the relative distance can either vary monotonically or can display 
oscillations, and the boundary between the two regimes is an unstable limit cycle. 

The implication of our results for the dynamics of biological organisms is twofold. First, we show 
that there are no stable equilibria (in position or orientation) between the cells, a result which is 
true arbitrarily of the sign of the far- field flow field each cell is generating (pushers or pullers). As a 



result, populations of ce 



Mendelson et al 



Sokolov et al 



2005 



1999 



2007 



Aranson et al 



Is are expected to a l ways 



Wu and Libchaber 



Cisneros et al 



2007 



2000; 



dynamically evolve, as is observed in experiments 



Dombrowski et al 



20071') and modeling 



Saintillan and Shelley, 



2007; 



2004 



Kim and Breuer 



Sim ha and Ramaswamv 



Ishikaw a and Pedlev 



2002 



2007b 



a; 



2004 



Hernandez-Ortiz et 



al. 



Saintilla n and Shellev 



2008 



Wolgemuth 



200 



Underhill et al 



2008 



Ishikawa and Pedlev 



2008; 



Mehandia and Nottl . 



20081 ) of cell populations, with an intermittence at the origin of the expression "bacterial tur- 
bulence". The model system studied in this paper allows us in particular to quantify rigorously 
the rate at which the cells are being effectively repelled from, or attracted to each other, and to 
obtain all types of possible swimming kinematics at t — > oo. In addition, what this paper shows, 
is that hydrodynamic interactions leads to "new" modes of swimming, meaning that the motion 
of each swimmer contains a component due to the presence of another cell which, over long times, 
integrates to an order one change in its swimming kinematics. This is reminiscent of recent work 



sho wing that hydrodynamic interac t ions can impart mot i 



ies ( Alexander and Yeomans 



2008; 



Lauga and Bartold . 



ity to otherwise non-swimming active bod- 



2008J), and is relevant to the experimental 
observation that dense cell p opulations display different length, time and veloc i ty scales than that of 



individual mic r o-org anisms (jMendelson et al 



Cisneros et al. 



20071) 



1999 



Dombrowski et al 



2004 



Sokolov et al 



2007 



5.2 Modeling assumptions and possible extensions 

The results in this paper were obtained under a number of simplifying assumptions, which we now 
discuss. 
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5.2.1 Non-spherical swimmers 



We have first assumed that the two swimmers are spherical, so that the rotation rate induced by 
the hydrodynamic interaction is equal to half the vorticity field created by the other swimmer. A 
corrective term of the same order however appears as soon as t he swimmer shape is not purely 



spher ical. Analytic solutions have been obtained for ellipsoids (jJefferv 



1991 



1922; 



Kim and Karillal . 



). In this paper, we focus on the spherical case as it provides the simplest system, and because 
it is a first good approximation of the shape of spherical organisms using cilia or flagella whose 
effect on the induced rotation rate can be neglected if their size is small compared to the body of 
the swimmer. If the organism is not spherical, a corrective term to the system, Eq. (|13p must be 
added to account for the effect of anisotropy and local strain rate. In the case of an e llipsoidal 



swimmer, the induced rotation rate on swimmer 2 is given by (jPedlev and Kesslei . 



1992) 



Ox^ 2 = i w W (r) + fa p 2 x (E« (r) • p 2 



(93) 



with r = r 2 — ri, and E 1 - 1 ) the vorticity field and strain rate tensor created by the motion 
of swimmer 1, p 2 the unit vector associated to the direction of the major axis of the ellipsoidal 
swimmer 2 and (3q = (c 2 — l)/(c 2 + 1) with c the ratio of major axis to minor axis of the ellipsoid, 
and measures the departure from the spherical case (p 2 moves rigidly with the swimmer). A 
reasonable approximation is to consider that p 2 ~ e 2 , i.e. the intrinsic swimming motion occurs 
in the direction of the major axis of the ellipsoid. The strain rate tensor is obtained from Eq. ([8]), 
and after substitution in Eq. (i93j) . the induced rotation rate becomes 



7 r x (S« -r) 



(94) 



+ 7^0 



5(r-S( 1 )-r)(e 2 -r)(e 2 xr) (e 2 • S« • r)(e 2 x r) (e 2 • r)(e 2 x (S« • r)) 



As pointed out above, the rotation rate now depends not only on the orientation of swimmer 1 
(through S^ 1 )) but also on the orientation of swimmer 2. In the limit of far-field interactions, the 
multiple-scale analysis of ^3.11 is still valid and we can study the average motion of the swimmers 
as represented by their mean distance r and the orientation of their intrinsic rotation vectors e^. 
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However, the relative phase between the circular motion of the swimmers (value between and 2ir) 
does not disappear in the averaged equations and acts as an additional arbitrary parameter. 



5.2.2 Validity of the far-field approximation and regularization 

When the two swimmers are not far enough from each other, two of our assumptions successively 
break down. Firstly, the separation of time scales is no longer valid when the time-scale associated 
with the intrinsic rotation of each swimmer is no longer much smaller than the time scale associated 
with the hydrodynamic interaction. In that case, the multiple scale analysis of $3] breaks down, 
and one needs instead to consider the full coupled equations, Eq. (fl"3j) . 

Secondly, when the swimmers get close to each other, the description of hydrodynamic interac- 
tions as being dominated by their far-field limit is no longer valid, and the following three terms 
need to be considered: (a) Higher-order corrections in the velocity and vorticity field created by a 
swimmer in Eq. (jlip (flows with r~ 3 decay such as force-quadrupoles, source-dipoles; flows with 
r -4 decay etc.). (b) Higher-order corrections in the induced velocity on a swimmer whose size is 
no longer negligible compared to the characteristic length-scale of the local flow. For a sphere, the 



( Jefferv . 


1922; 


Lamb. 


1932) 



1932 ). For arbitrary shape, ge neral frameworks have been studied allo wing 



1964 



Liron and Bartal . 



1992J); (c) 



the computation of the successive corrective terms (jBrennerl . 
Higher-order corrections due to the two-way coupling: Swimmer 1 creates a flow field that influ- 
ences swimmer 2. But the presence of swimmer 2, modifies this flow field (even if swimmer 2 was 
not swimming) which also induces a correction on the velocity of swimmer 1. 

These three contributions are negligible for large distances but can become dominant at in- 
termediate distance or in near-field interactions. In particular, we observed previously that the 
far-field behavior can lead to collisions between the swimmers (when a or r go to zero) as the 
hydrodynamic interaction terms in Eq. (|13|) are attractive for particular relative orientations of the 
swimmers, regardless of their relative distance. Moreover, these attractive interactions also have a 
diverging amplitude as r — > 0. Obviously, the far-field approximation is violated when the distance 
becomes small, and the higher order corrections discussed above must be included to account for 



regularizing forces that arise at intermediate or short distances. In an effort to remain general, one 
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could attempt to introduce some empirical short distance corrections to reduce the attractive terms 
in the near- field (in the form of exponential or power-laws regularization at small r for example), 
but these are not based on physical principles. For distances between the swimmers much smaller 
than the typical size of each swimmer, lubricati on theory can be use d but, in the intermediate 



distance range numerical simulation is necessary (jlshikawa et al 



2006) 



In general, if one is interested in intermediate or short range interactions, a knowledge of 
the detaile d swimmer geometry and propulsion method is n ecessary, as is the case for spherical 



squirmers (llshikawa et al 



( Hernandez-Ortiz et al 



200 



2006; 



Ishikawa a 



Gvrva et al 



nd Pedlevl. 



2007b 



a|) or dumbbell-like model organisms 



20091) • Squirmers maintain a spherical shape at all time 



and generate motion by tangential displacement of their surfaces. They are generally thought as 
a good approximation for spherical swimmer using ciliary propulsive schemes, the spherical shape 
of the swimmer corresponding to the envelope of the cilia in that case. The squirmer formulation 
has the advantage that an analytic solution exists for the velocity field created. Using Faxen's 
law, an exact system for spherical swimmers can then be obtained. Analytic solution of the multi- 
body problem is however not possible in general and such a system must be solved numerically. 
Both squirmers and dumbbell-like organisms are simple approximations of real swimmers, and 
considerations of the detailed geometry of the swimmer often lead to a trade-off between accuracy 
in the biophysical description of real organisms and simplified representations to allow an easier 
mathematical or numerical treatment. 
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